Modifications for numerical stability of black hole evolution 
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We experiment with several new modifications to the Baumgarte-Shapiro-Shibata-Nakamura 
(BSSN) formulation of Einstein's field equation, and demonstrate how these modifications affect 
the stability of numerical black hole evolution. With these modifications, we obtain accurate and 
stable simulations of both single excised Kerr-Schild black holes and punctured binary black holes. 
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I. INTRODUCTION 

Numerical relativity is aimed at solving Einstein's 
equations with the aid of computers. It took decades to 
reformulate Einstein's equations to the required stabil- 
ity and accuracy in simulations. Breakthroughs in 2005 
and 2006 |1, 2] brought this development to a more ma- 
ture status, gave confidence to the community in mod- 
eling gravitational wave sources and extracting informa- 
tion of gravitation radiation from the simulations of bi- 
nary compact objects. Since then, observations in simu- 
lations which have been refined and extensively studied 
over the past few years include long term gravitational 
waves and the final state of the binary compact objects 
merger (see review articles ■'> "j and reference therein), 
gravitational recoil of a black hole m, relativistic jet for- 
mation from mergers of black holes [7| , and neutron stars 
Q. These also provide new insights into mathematical 
general relativity. Numerical relativity has now become 
an indispensable and efficient tool in the research of gen- 
eral relativity and astrophysics. 

Among the numerous reformulations of Einstein's 
equations, two particular formulations are most fre- 
quently adopted for the simulations of black holes. One 
of them is the generalized harmonic formulation, in cither 
second-order formulation Q, or fully first -order formula- 
tion [9( (In either case the key ingredient for stability of 
this formulation is a constraint damping mechanism 10]). 
The other is the Baumgarte-Shapiro-Shibata-Nakamura 
(BSSN) system [ll|, which has been implemented by 
many groups using finite difference codes, in first-order- 
in-time and second-order-in-space form. 

There have been many studies in modifying the BSSN 
formulation to increase its numerical stability. For exam- 
ple, Alcubierre and Brugmann [l2[ combined the BSSN 
formulation with constraints enforcing the traceless con- 
dition of the conformal extrinsic curvature in every time 
step, replacing the conformal connection function T l with 
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that calculated from the Christoffcl symbols whenever 
r* is undifferentiated. With these modifications, they 
evolved single black holes in Kerr-Schild coordinates sta- 
bly with octant grid symmetry, and all fields settled down 
to equilibrium without encountering any instability. The 
confinement to grid symmetry was relaxed later in [ilij . 
by adding the T-constraint to the field equation of P 
in order to suppress instability, as well as by employing 
alternative techniques in the enforcement of other alge- 
braic constraints. Yoneda and Shinkai analyzed analyti- 
cally the constraint propagation of BSSN formulation in 
[lij , and proposed adjustments to obtain better stability 
by adding constraint terms. The follow-up work in [151 ] 
numerically verified the advantage of their adjustments 
in stability over the original BSSN formulation. Higher- 
order derivatives of constraint was also added to the field 
equation in [l6| to enhance stability. A first-order BSSN 
formulation has been developed in [l7| to seek more sta- 
ble performance in simulations. 

In this work, we report numerical tests of new modifi- 
cations to the BSSN system. The work can be considered 
as an extended study based on an earlier article . The 
modifications include (1) adding the T-constraint to the 
field equation of the conformal three-metric; (2) replac- 
ing the conformal connection function calculated from 
the conformal metric, i.e., r g , with the independent T l , 
and enhancing the derivative of the unimodular deter- 
minant constraint in the connection with an irreducible 
decomposition; (3) applying, with some deformation, the 
adjustment of the field equation of the conformal extrin- 
sic curvature with the momentum constraint proposed in 
We also emphasize the alternative method in [l3[ 
in the enforcement of the unimodular constraint and the 
traceless conformal extrinsic curvature constraint in ev- 
ery timestep. We experiment with these modifications on 
simulations of single Kerr-Schild black hole, and obtain 
evolutions with long-term stability. These modifications, 
also applied to the evolutions of binary punctured black 
hole, demonstrate better stability and accuracy than the 
original BSSN formulation against numerical errors, ei- 
ther from finite-differencing or from mesh refinement. 

The paper is organized as follows: We summarize the 
BSSN formulation in Sec. II. Our modifications of the 
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BSSN scheme are described in Sec. III. Single black hole 
spacetimes in Kerr-Schild coordinates and binary black 
hole with punctures are described in Sec. IV. We discuss 
numerical implementations and the gauge conditions in 
Sec. V. In Sec. VI we present results of our simulations 
for both single black hole and binary black holes. We 
summarize and discuss the implications of our findings 
in Sec. VII. Throughout this paper we adopt geometric 
units, G = c = 1. 



with the conformal exponent <f> chosen so that the deter- 
minant 7 of 7y is unity 



(10) 



where 7 is the determinant of 7^ . The traceless part of 
the extrinsic curvature K^, defined by 



Aij = K {l]) = Kij - -j tJ K, 



(11) 



II. THE BSSN FORMULATION 

The metric in the ADM form is 

ds 2 = -a 2 dt 2 + j l3 {dx l + p&t^&x 3 + (5h\t), (1) 

wherein a is the lapse function, f3 l is the shift vector, 
and jij is the spatial three-metric. Throughout this pa- 
per, Latin indices are spatial indices and run from 1 to 
3, whereas Greek indices are space-time indices and run 
from to 3. 

Einstein's equations can then be decomposed into 
the Hamiltonian constraint H and the momentum con- 
straints M.,, 



H = R — K i: jK ij + K 2 = 0, 
i 

and the evolution equations 



M, = VjK'i ~ ViK = 0, 



— 7y = -2aK v 



(2) 
(3) 



(4) 



—Kij = -ViVjCt + a(Rij - 2K u K l t + KK tj ). (5) 
Here we have assumed vacuum T Q « = and have used 



P 

dt dt P 1 



(6) 



where £^ is the Lie derivative with respect to (3 l . V, 
is the covariant derivative associated with 7^ , R^ is the 
three-dimensional Ricci tensor 



1 



Rij ~~2^ (^ k ^ a 1fit,kj — lH,ij — lij,ki) 



+ 7 (r m «r 



where 



(7) 



(8) 



And R is its trace R = Y^Rij ■ 

In the BSSN formalism [ll| , the above ADM equations 
are rewritten by introducing the conformally related met- 
ric 7ij 



(9) 



where with two indices between () is to take the sym- 
metric and traceless part of Kij, and K — j 1 ^ Kij is the 
trace of the extrinsic curvature, is conformally decom- 
posed according to 



Aij — 6 ^ Aij . 



(12) 



The conformal connection functions T l , initially defined 
as 



(13) 



are regarded as independent variables in this formulation. 

The evolution equations of BSSN formulation can be 
written as 



dV = -^' 






(14) 


d - 2 
— 7y = -2aAij, 






(15) 


±K = a(AijA* 




- V 2 a, 


(16) 


—Aij = a(KA l3 - 


2A lk A k j) - 


h e-^(aR m 


-v (l v j)a ), 






(17) 


d t V = 2a (f)^ 


3 ' 




- 2A ij aj 


+/3 j f\j - f j /3 l 






(18) 



The Ricci tensor Rij can be written as a sum of two 
pieces 



Rij = Rij + R^, 



(19) 



where Rfj is given by 



Rfj = -2ViV^ - 27yV 2 (/. + 4Vi«£Vj<£ - ^ijV k 4>V k 4>, 

(20) 

where V; is the covariant derivative with respect to 7^, 
while, with the help of the T l , Rij can be expressed as 



Rij = - -ji mn iij,mn + 7fe(»r j) + r fc r (iJ -) fc 

+ 7 mn [2r m {iXj)kn + F inFkmj]. (21) 
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The new variables are tensor densities, so that their Lie 
derivatives are 



£ K 


= P k K tk , 




(22) 




= P k 4>,k ■+ 




(23) 




= P%j,k 


+ ,j) - -^lij(3 k ,k, 


(24) 


£ pAij 


= f3 k Aij^ 


+ 2^4fe(i^ fe j) - -A l j(3 k yk . 


(25) 



The Hamiltonian and momentum constraints © and ([3]) 
can be rewritten as 



H = e-^iR - 8V 2 </> - 8VfyVi0) + -K 2 - = 0, 

(26) 

(27) 



M = VjAS + 60, AS - -K i = 0, 



where R — j v Rij. Besides being used to obtain the evo- 
lution equations (H"6|) and (H"8l) in the BSSN formulation, 
the Hamiltonian and the momentum constraints are also 
applied to the volume integrals of the ADM mass and 
the angular momentum (in vacuum), respectively [l8j : 



M = 777" i (r 1 ~ m^dje^dti 
lb7T Jan 



l - J lefV&jA* - \k 2 ) + - e*R]d 3 x, 



16tt Jn 
1 

8^ 
1 

'87^ 



8tt Jan 



(28) 

(29) 
(30) 



(i^fe + -x>K ik - -x^A em ^ m >k )d 3 x, 

(31) 



where dS^ = (l/2)eyfcda; J dar\ These two global quanti- 
ties are useful tools for the system diagnostics to validate 
the calculations. 



III. ADJUSTING THE BSSN EQUATIONS 

The BSSN formulation has been described in detail in 
previous papers We will discuss here only the new 
improvements. For a solution of the BSSN equations to 
be equivalent with a solution of the ADM equations, the 
new auxiliary variables have to satisfy new constraint 
equations. In particular, the determinant of the confor- 
mally related metric 7^ has to be unity, 



D = 7-l = 0, 



and An has to be traceless 



T = f J A l3 =0, 



(32) 



(33) 



and the conformal connection functions T 1 have to satisfy 
the identity 



ff i = r i -r i =f i + f i j = o l 



where T l „ = ^ k T l ik- These conditions 



jk 



(34) 

(|34l) are also 



viewed as constraints in the BSSN formulation, in addi- 
tion to the Hamiltonian and momentum constraints. It 
is worth mentioning that in the recent efforts on the for- 
malism extending the solution space of Einstein's equa- 
tion the G l 's are related to new dynamical vari- 
ables whose evolution is mainly driven by momentum 
constraint, and therefore can be regarded as the cumu- 
lated effect of momentum constraint violations. 

In an unconstrained evolution calculation, the con- 
straints are monitored only as a code check. It has been 
proven to be advantageous, however, either to enforce at 
least some of the constraints during the evolution, or to 
add evolution constraint equations to the evolution equa- 
tions. 



A. Enforcement of the constraints T> and T 

In the conventional adjustments [H|>[2(|, the algebraic 
constraints (|32|) and (j33|) are enforced actively by replac- 
ing jij and Aij with the following: 



A, 



A 



(ij) ' 



(35) 



after every time step. With Eq. (1351) . the noise from the 
violation of constraints (|3^|) and ([531 is scaled/subtracted 
"evenly" from each conformal metric/extrinsic curvature 
component. These two adjustments are widely used in 
most of the numerical relativity groups. 

The alternative adjustments for constraints (|32p and 
([33)1 in are as follows: Instead of treating all compo- 
nents of 7y equally, only five of the six components of 7y 
need to be evolved dynamically, and the remaining one 
can simply be computed using Eq. (|3"2"]) . For example, let 
"f zz be the chosen component. Then 



' lyylxz 



^Ixylyzlxz + Ixx 7 



^/xx^fyy 



>xy 



(36) 



where % x , ^ yy , % y , j yz , j xz are evolved with the field 
equation (I15[) . In principle, any one of these six compo- 
nents of jij can be chosen to be computed using Eq. (|32l) , 
leaving the other five to be evolved with Eq. (fl5|) . How- 
ever, there will be extra difficulty if any of the three 
off-diagonal variables is chosen since Eq. (|32|) gives a 
quadratic equation for an off-diagonal component instead 
of a linear equation for a diagonal component. Similarly, 
only five of the six components of A^ need to be evolved 
dynamically, and the remaining one can simply be com- 
puted using Eq. (|33)) . For example, let A yy be the chosen 
component. Then 



A): 



Ax* + A Z Z + ~A X yl XV + Ayzl V% 



(37) 
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Although any one of these six components can be cho- 
sen to be computed using Eq. (1551) , leaving the other five 
to be evolved with Eq. (fT7| . However, it is not recom- 
mended to choose any of the three off-diagonal compo- 
nents since the corresponding denominators for the off- 
diagonal components, i.e., j xy , 7 yz , and 7 XZ , could vanish 
anywhere. 

One of the features of the alternative adjustments on 
jij and Aij is the economy compared with conventional 
methods. Only five, instead of six, components for both 
jij and need to be evolved dynamically, and the re- 
maining one is determined by the algebraic constraint. 
Meanwhile, the alternative adjustments "correct" only 
one component instead of all the components. It can be 
expected that results with the alternative adjustments 
will be more convergent than those with conventional ad- 
justments. On the other hand, the obvious shortcoming 
for the alternative adjustments is the asymmetric treat- 
ment of the six components. However, as far as the cases 
we have ever checked, the effect of the asymmetry is neg- 
ligible. It is helpful numerically to choose a different 
diagonal pair for 7^ and Aij , instead of the same diago- 
nal pair, in eqns (|36[) and (|37p to increase the asymmetry 
in the evolution equations and to suppress the possible 
growth of the unstable modes ignited by numerical er- 
ror. In the work we choose the pair (7 ZZ , A vv ), instead of 
(7z2, A~zz) in [H| , to fulfill this requirement. 

B. Decomposition 

For a third-rank tensor Yijk with symmetry in the last 
two indices, i.e., Yak = Yi(jk)i it can be decomposed into 
the following form [21j 

3 11 

Y tjk = ¥ ijk + -li(jL k) - -li(jY k) + -jjkYi, (38) 

where 

Li = ~f jk Y jki = j jk Y jik , Y = 7^'% fe , (39) 

and ¥ijk is the traceless part of Yij k , i.e., ¥ik h = ¥ k ik = 
¥ fc fei = 0. We can apply the decomposition (|38l) to the 
connection in every time slice although a connection is 
not a tensor. The is because any quantity can be decom- 
posed like a tensor as long as the quantity is only con- 
sidered in the same coordinate, without any coordinate 
transformations involved. Thus the conformal connec- 
tion can be decomposed as 

f*i* = Pik + \#<jf h) - ^ (i ffc> + \lA ( 4 °) 

where Ti = T k ki — ^ k ik, and F % j k is the traceless part of 
r jk . Here % = T k kl = djln V7 = analytically, but Tj 
could be nonzero numerically due to the truncation error. 
We can guarantee the vanishing of T$ by subtracting the 
terms having it from Eq. (14T))) . T l in Eq. (|4T))) can be 



replaced with the conformal connection function f * by 
adding the T-constraint to it. Then the new connection 
becomes 

=**ik - l^ijh) + p jk f \ (41) 

We substitute T z j k with the new one T % j k in all the cal- 
culations, including the conformal covariant derivatives 
and the conformal Ricci tensor, as a new modification. 

Similarly, the spatial derivative of the conformal metric 
can also be decomposed in every time slice as 

3-1 1 

dij jk = da jk + -7*0-1% - g7i(7 r *) + -ftjkTi, (42) 

where dijj k is the traceless part of 9j7jfe. We use the 
same trick on the conformal connection part like we did 
on Eq. (|4ip and obtain a new spatial derivative on 7^ 

3 

d'iljk = dtfjk + -^li(jGk) ■ (43) 

Here we do not take action on eliminating the T l part in 
Eq. (|4"2"| since its effect is negligible from the observation 
of our numerical experiment . 

However, it turns out that this replacement of the spa- 
tial derivative of the conformal metric seems give too 
much change on the field equation of 7^ and causes in in- 
stability when d'ijjk is applied in Eq. ([Ml) . This problem 
can be solved when the modification ([56]) in Sec. IIIIDl is 
applied in simulations, at least in single black hole simu- 
lations. Nevertheless, we modify Eq. (|43p . by multiplying 
an adjustable parameter to the substitution part, to have 

dajk = dijjk + &7i(jGk) - \ljkQi, (44) 

where the range of the parameter is usually chosen as 
a e [1/5,4/5] in this work. And 

Fdajk = Pdajk + crPuQk) - pjkP%- (45) 

Thus the replacement of % di*jj k in Eq. (l2~4l with /3 l di^jk 
is equivalent to modifying Eq. (1151) into 

^7y = -2aA ij + a/3 (i g j} - pi^ k Q k . (46) 

We postpone the modification on the field equation of 
Aij with the decomposition of the spatial derivative of 
Aij until Sec. IHlD1 

C. Enforcement of the T-constraint 

In the conventional adjustments, to enforce constraint 
(1M1) all the undifferentiated T 1 in the evolution equations 
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are substituted with T l g . However, this adjustment give 
stability in single Kerr-Schild black hole simulations only 
in octant symmetry [l2| • Instead of the conventional ad- 
justment, one of the alternative adjustments [l3| (dubbed 
as "YBS") is to add the T-constraint to the evolution 
equation JTS]) of F by 

d t r = rhs of (HHJ)- |(e + l)^/5 fc ,fe, (47) 

where £ is usually chosen to be unity. The YBS adjust- 
ment has been proven to be helpful in suppressing the 
instability caused from some unstable modes [HI, [22| ■ 

Here we propose another adjustment to enhance the 
stability. This adjustment is basically the following sub- 
stitution in Eq. (fT8|) : 

f j Kj -> K ij j +KT\ (48) 

where k 1 ^ = j^K. And this adjustment turns Eq. dT5J) 
into 

d t f 1 =2a(f i ife A fe - ^ K ij d + 6i w 0j) - 2A ij a tj 

+ ftf\j - t^\j + ^ k p\ jk + p ij p k , jk 

+ ^(/3 fe , fc - 2aK)f\ (49) 

As we will show in Sec. I VII the performance of the code 
with this adjustment is better than the one with the YBS 
adjustment, but without the uncertainty of choosing the 
value of the parameter £. Furthermore, the connection in 
Eq. (j49|) can be replaced with the new connection r % j k , 
and thus Eq. (|4"9"]l is modified into 

d t f l =2a(r i jk A jk - \k 1 \j + GA^fa) - 2A^a d 

+ ^r.j - + i jk p\ jk + p ij P k , jk 

+ ^{fi\ k - 2aK)r - (1 + £)@(X i )X i G i , (50) 

where the last term is newly added to control the stability 
via the linear term of T* , is a step function, and \ l is 
as follows: 

\ l = \{P\ k ~2aK)-P ;i -^ a Al (51) 

where the index with hat, i.e., i, means that no index 
summation happens on this index. £ is chosen to be 1 in 
all the cases in this work. 



D. Application of the momentum constraint 

Yoneda & Shinkai have studied the adjusted systems 
for the BSSN formulation in It shown in the work 
that the adjusted BSSN system could be quite robust 



with the following modification on the field equation (fl7|) 
of Ay : 

-^■Ay =rhs of (PU) +K A aV(iMj), (52) 

where ka is a constant. If ka is set as positive, the vi- 
olations of the constraints are expected to be damped. 
The robustness of this Aadjusted BSSN formulation has 
been demonstrated in [l5|, [Tg . Here we would like to ap- 
ply this type of modification to the BSSN formulation in 
a slightly different manner. In the momentum constraint, 
the covariant derivative of Ay can be rewritten as 

VjAV = Aj J t j — r ^ i-A k j = Ai\j — jk.i- (53) 

Thus the momentum constraint turns to be 

Mi = Ajj - \^ k A jKl + 64> tj A^ - ^K ;i . (54) 

And the symmetric part of its spatially partial derivative 
is 

=A(i h ,j)k + f U (i\Aki,\j) - -jl M A k t,ij 

~ t - i 2 

+ Q&XiAjf + 6<^A(i ,j) ~ l K ,ir (55) 

Therefore, the modification on Eq. (IT71) in this work is 
—Ay = rhs of (HZD + hf(a)M {itj) , (56) 

where h is the grid width and f(a) is a function of lapse 
and usually chosen to be 1 in single black hole simu- 
lations. It shows in our numerical experiments that this 
modification is helpful in suppressing the instability from 
the high-frequency unstable modes. 

As in Sec. MI B[ the spatial derivative of Ay can be 
decomposed in every time slice as 

— - 3 11 

8iA jk = diA jk + -li(jP k ) ~ g7i<jQfc) + ^IjkQi, (57) 

where diAj k is the traceless part of diAj k and 

Pi = y k diA ki , Q, = f k diA jk . (58) 

With the momentum constraint (I5~4l and the traceless ex- 
trinsic curvature constraint (1331) , we obtain a new spatial 
derivative on Ay 

3 1 1 

diAjfc = diAjk - -ji(jM k ) - y^7ioAfc} - ipjkAi, (59) 

where A\ is the spatially secondary constraint of Eq. (|33[) 

A = diT = 7 jk A jkii + A jk f k tl = Qi- 2A jk r^ k i = 0. 

(60) 

With the new spatial derivative (f51?|). we can modify the 
field equation of Ay further from Eq. (fB"6"|) to 

d 3 
—Ay = rhs of (HZD + hf(a)M {iJ) - -(3 {l M 3) 

- ^a A i) - pi^Ak. (6i) 
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IV. INITIAL DATA 

A. Single black hole in Kerr-Schild coordinates 

The ingoing; Kerr-Schild form of the Kerr metric is 
given by [23ll24| 



ds 2 = (r)^ + 2H£^ u )Ax' x &x v , 



(62) 



where rj^ u — diag(— 1, 1, 1, 1) is the Minkowski metric 
in Cartesian coordinates, and H a scalar function. The 
vector is null both with respect to r]^ and g^ Ul 



In the static case a = 0, the above expressions reduce 
to the Schwarzschild expressions in ingoing Eddington- 
Finkelstein form [2a| 



r 



2M 



(1,^), 
r 



Kij r 4 (l + 2M/r)V2 



K = 



2M 



r 2 (l + 2M/r) 3 / : 



r r]ij - (2 + —)xiXj 
■(1 + 3—), 



(73) 



where M is the total mass-energy and r 2 = x 2 + y 2 + z 2 



(63) 



B. Binary black hole with punctures 



and we have £ 2 = l k l k . The general Kerr-Schild black 
hole metric has 



H = 



Mr 



• 2 + a 2 cos 2 



and 



rx + ay ry — ax z 
r 2 + a 2 ' r 2 + a 2 ' r 



(64) 



(65) 



Here M is the mass of the Kerr black hole, a — J/M 
is the specific angular momentum of the black hole, and 
r and 9 are auxiliary spheroidal coordinates defined in 
terms of the Cartesian coordinates by 



V 



(66) 



and z = rcosO. The event horizon of the black hole is 
located at 



M + a/M 2 - a 2 



(67) 



Comparing (|62|) with the ADM metric (JTJ one identifies 
the lapse function a, shift vector ft and the spatial 3- 
metric 7y as 



1 



ft = 2#4 
7u = Vij + 2H£i£j. 



(68) 

(69) 
(70) 



We can see here that these variables all extend smoothly 
through the horizon and their gradients near the hori- 
zon are well-behaved. Given these metric quantities, the 
extrinsic curvature can be computed from (£[]) 

K l3 = 2aHt k [i i £jH ik + 2He {i £ jhk ] 

+ 2a[£ {i H j) +m (iij) }, (71) 
K = 2a\l + H)l k H M + 2aH£ k , k . (72) 



We use the quasi-circular binary black hole puncture 
data to test our modifications. The momentum param- 
eter for the quasi-circular orbit is set to be the value 
given by (2(| based on the helical Killing vector condi- 
tions. In the following, we review the puncture scheme 
and describe how we construct initial data by the multi- 
domai n sp ectral method extended from the LORENE li- 
brary Ha S|. 

To determine a three-geometry subject to the con- 
straint equations and ([3]), the conformal decompo- 
sition plays an important role (see [29j j for instance). 
Lichnerowicz [30l | proposed a conformal decomposition 
of three-metric 7^ = ip 4 ^, as Eq. with ip = e^, and 
later York [3l[ used the transverse-traceless decomposi- 
tion of the conformal extrinsic curvature A 1 - 7 = tp^A 1 ^ = 
A l ^ T + LW lJ . The particular conformal scaling makes 
the identity ; V i A ij = -0" 1O V \A 11 hold. In this decom- 
position, A l ^ T is transverse (i.e., ViA l ^ T = 0) and can 
be assumed to be zero to have a purely longitudinal 
A ij = W' 3 , In this approach, named as conformal 
transverse-traceless method (CTT), and in the assump- 
tion of conformal flatness, 7^- = rjij, maximal slicing, 
K = 0, and A^ T = 0, the vacuum constraint equation 
would then be expressed as 



1 



7 A- ■ A l 3 



0. 



(74) 
(75) 



The momentum constraint is linear in this case. For a 
iV-black hole system, it allows the Bowen-York solution 
j^-j = Yl a =i a j where the conformal extrinsic curvature 
for each black hole is 1321 



A% = ip 10 A*i = — -[pCW* - 2(j ij - n l n j )P k n k ] 
At 1 

+ ^n^e^ M S k n e , (76) 
2r A 

with n % the spatial unit vector pointing away from the 
puncture, and P l and S l its linear and intrinsic angular 
momentum respectively. The conformal factor then can 
be numerically solved from the Hamiltonian constraint. 
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TABLE I: Input parameters and modifications used for selected evolutions. For each evolution we list the symmetry used, 
the method for constraints T> & T enhancement, the modifications on the field equations of Y l and 7^, the substitution of 
connection, the modification on the field equations of Aij, the time when instability shows. 



Case 


Symmetry 


Jzz+Ayy 


r 






Aij 


Instability 


Eqn 






(EH 


(Snj) 


<EU 


(EH 


appears at 


STD 


Equ 


([35]) 










~ 250M 


YBS 


Equ 


V 


(S3 








~ 650M 


El 


Equ 


V 










~ 550M 


E2 


Equ 


V 




cr=2/5 






~ 750M 


E3 


Equ 


V 


V 


a=2/5 






~ 850M 


E4 


Equ 


V 


V 


a=2/5 


V 




~ 950M 


E5 


Equ 


V 


V 


a=3/5 


V 


(El 


None 


E6-I0 


Equ 


V 


V 


a=3/5 


V 


V 


None 


E6 


Equ 


V 


V 


a=3/5 


V 


V 


None 


E6-hi 


Equ 


V 


V 


a=4/5 


V 


V 


None 


NT 


None 


V 


V 


a=3/5 


V 


V 


None 



In the puncture method, one separates out the singular 
part 

N 

^ = 1 + E^ (^) 

a—l 

from the conformal factor ip, where m a and r a are the 
mass parameter and coordinate distance from each punc- 
ture respectively. It is the superposition of Schwarzschild 
solution in the isotropic coordinate and satisfies the 
Laplace equation [33| . Therefore the desired regular part 
u = 1/) — tp s is determined by 

V 2 u^~\A^A l3 ^ s + u)- 7 . (78) 
o 

The existence and uniqueness of the solution has been 
discussed in [33]. 

Our multi-puncture initial data solver is motivated 
from [35| for the excised binary black hole initial data. 
We cover on each black hole a spherical multi-shell do- 
main, and split u — u a (the index a runs over the 
number of the black holes) as well as the puncture equa- 
tion (J7HJ) into 

V 2 u a = -^A tJ (ij s + u)- 7 . (79) 

Only one of the extrinsic curvature tensor Aij split in the 
above equation. Thus it is expected that the source term 
in its RHS has large contribution only near each hole, and 
the use of spherical polar coordinates is adequate for solv- 
ing the equation near the punctures. The u a vanished at 
the asymptotically flat physical outer boundary on the 
outermost, compactificd shell. And d r u a (r = 0)= is 
also ensured at the punctures as the inner boundary con- 
dition. These equations for each hole are then solved it- 
eratively with the Poisson solver in the LORENE library 



until each successive difference of Su a is small, typically 
1CP 11 . For the binary case, our initial data is consis- 
tent with the earlier wor k 13611 . The convergence of our 
method was presented in [371 ]. 



V. NUMERICAL IMPLEMENTATION 

In this work, we discretize the evolution equations us- 
ing a fourth-order Runge-Kutta scheme. We use fourth- 
order centered differencing everywhere except for the ad- 
vection terms on the shift. For these terms, a fourth- 
order upwind scheme is used along the shift direction. 

For the lapse gauge condition, we consider the "1+log" 
slicing [38l440l ] which is basically a modification of the 
Bona-Masso slicing. The lapse condition used in the sin- 
gle Kerr-Schild black hole evolutions is 

d t a = D t f3 l - aK. (80) 

For binary black hole simulations, the moving puncture 
technique is adopted. And the lapse condition used here 
is the standard one in numerical relativity community for 
binary black hole simulation as 

dta = 0*a,i - 2aK. (81) 

Many driver g aug e conditions (e.g., the T-driver) for 
the shift vector [4ll. l42j are currently the main type of 
gauge conditions used in the punctured black hole calcu- 
lations. In this work, we will only focus on these types of 
gauge conditions. The hyperbolic type T-driver condition 
used for the single black hole simulations is 

dtP 1 = \b 1 , (82) 
d t B l = a 2 d t f l -7]B\ (83) 
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where r/ is the parameters to be chosen. In the single 
Kerr-Schild black hole evolutions, r\ = 5. For the binary 
black hole simulations, we use 

8 t i = lB i + p>'l3 i j, (84) 
d t B* = dtf* - 2B l + p j B\ j - /3 j f\j. (85) 

On the outer boundaries of the numerical grid we im- 
pose a radiative boundary condition that is imposed on 
the difference between a given variable and its analytic 
value / — /analytic = u(r — t) / V where u is an outgoing 
wave function. We apply this condition to all fields ex- 
cept P which we leave fixed to their analytic values at 
the boundary in the single Kerr-Schild black hole evo- 
lutions for stability. We return to apply the radiative 
boundary condition on T l in the punctured binary black 
hole evolutions. 

The excision technique is applied in the single Kerr- 
Schild black hole evolutions since the singularities in the 
Kerr-Schild coordinates are physical ones for which the 
puncture method is not applicable. (Usually the punc- 
ture method is applied to the isotropic-like coordinates in 
which the singularities are coordinate ones.) For (spher- 
ical) excision regions, we adopt the recipe suggested by 
[l2l Il3j to copy the time derivative of every field at the 
boundary from its neighboring grid-point. The details 
can be found in Hal. 




, 1 , , . , , , , , E? , , , , , ,E5, , 
500 1000 1500 2000 

t/M 



FIG. 1: The root mean square (rms) of the change in the 
trace of extrinsic curvature between consecutive time steps as 
functions of time for the cases listed in Table [T] The usual 
BSSN formulation is used in Case STD. The modifications 
suggested in [l3|] are used in Case YBS. The modifications 
suggested in this work are activated one by one from Case 
El to Case E6. Case E5 has the same modifications with 
Case E4 except with a different value of a used in Eq. (|46|) . 
Case N7 has the same modifications with Case E6 except with 
none symmetry for the computational domain. There is no 
instability shown in Cases E5, E6, and N7 throughout these 
runs. 



VI. NUMERICAL RESULT 

A. Single black hole 

Our simulations for single Kerr-Schild black holes are 
summarized in Table HI For most of these simulations we 
use computational domains of size — \2M < x,y < 12M 
and < z < 12M for equatorial symmetry, and — 12M < 
x,y,z < 12M for no symmetry, with a grid spacing of 
Ax = 0.2M. In order to analyze the effect of resolution 
we performed the two Cases E6-lo and E6-hi on the same 
domain and used a resolution of Aa; = 0.4M for E6-I0 and 
Ax = 0.1M for E6-hi. We always use a Courant factor 
of 1/4 so that At — A.x/4. We excise spheres of radius 
r = 1.6M inside the event horizons (r = 2M) in these 
simulations. 

A simulation can be judged as being stable if changes 
in all dynamical variables drop to round-off error (of 
about 10 -16 for double precision), and remain at that 
level as the simulation goes. Reaching round-off implies 
that the numerical solution has settled down to the equi- 
librium solution of the finite-difference equations (as op- 
posed to the equilibrium solution of the differential equa- 
tions, which is provided as initial data). In all our stable 
runs, besides monitoring the global quantities consisting 
of the ADM mass ([251 |2"5|) . the angular momentum J z 
(|3QI3ip . the L2 norms of the Hamiltonian constraint %, 
the momentum constraint J\4 X , and the T-constraint Q x 



violations, we also monitor the changes in the represen- 
tative variables, i.e., <p, K, % x , A xx , T x , a, and (3 X , to 
see if they can reach and remain at the level of round- 
off error. Due to the existence of singularity inside the 
excision area, we only consider the domain outside the 
excision (including the part inside the horizon) for the 
estimate of the error norm and the changes. 

In Fig. [T] we summarize the listed cases with the 
medium resolution in Table HI As shown in [l2l . [l3| . 
in Case STD the evolution with the usual standard 
recipes applied to the BSSN formulation becomes unsta- 
ble quickly. The root mean square (rms) of the changes 
in K between consecutive time steps, A/f rms , in this case 
drops exponentially until t ~ 250M, then at later times 
it increases exponentially. This exponentially growing 
mode can be extrapolated back to about round-off er- 
ror at t = 0, indicating that the mode is triggered by 
round-off error in the initial data. In Case YBS the 
recipes suggested in [l3j], including mainly Eqs. (|36I37I) . 
and Eq. (T4"T|) . are applied to the BSSN formulation. Un- 
doubtedly, the evolution in Case YBS is more stable than 
the one in Case STD. We can see it from that the A-K" rms 
in this case drops exponentially until t ~ 650M before it 
turns to increase exponentially. In fact, the recipes used 
in Case YBS have been shown in 13] to be able to stabi- 
lize an evolution until its AA" lms reaches round-off error. 
However, a lower grid resolution and the second-order fi- 
nite differencing methods both spatially and temporally 
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FIG. 2: The monitored quantities as functions of time for Case E6. The upper-left panel compares different integrals for the 
ADM mass. The lower-left panel compares different integrals for the angular momentum. The upper-right panel shows the L2 
norms of the Hamiltonian constraint H and the momentum constraint M x . The lower-right panel shows a log plot of the rms 
of the changes in the lapse and the trace of extrinsic curvature between consecutive time steps. 



are used in those numerical experiments. The instability 
shown in Case YBS indicates that higher-frequency un- 
stable modes triggered with higher grid resolution and/or 
higher-order finite differencing methods need to be tamed 
with further modifications to the BSSN formulation. 

We test our new modifications from Case El to Case 
E6. In Case El, the field Eq. (gSJ for f\ instead of 
Eq. (|T7|) in Case YBS, is employed in the evolution. As 
described in Sec. IIII C| this modification is to enhance 
the effect of the linear terms in the field equation of P 
on convergence and stability. We can see from the plot 
that the evolution in Case El converges faster than the 
one in Case YBS, and thus encounters the growth of in- 
stability from round-off error earlier (t ~ 550M). Besides 
the modifications employed in Case El, we add the mod- 
ified field Eq. (|46p of 7^ in Case E2 with the parameter 
choice a = 2/5. In our numerical experience, this modi- 
fication is a critical one among the new modifications on 



the stability of an evolution. The choice of the parame- 
ter a could also affect the behavior of convergence. From 
Fig. [JJ we can see that the evolution in Case E2 does not 
show any instability until t ~ 750M. In Case E3, the 
connection T l jk is substituted by the new one, r l jk, via 
the application of Eq. (l4Tj) . The purpose of the modifi- 
cation is to substitute T l g , a decomposition part of T l jk, 
in F l jk with the independent conformal connection T l , as 
well as enhancing the secondary constraint T^j ~ 0. We 
can see in Fig. [T] that the combination of the modifica- 
tions used in Case E3 enhances quite a lot the stability 
of an evolution. The evolution does not show any in- 
stability until t ~ 850M and the A_ftT rms reaches almost 
round-off error before the unstable mode appears. 

However, we found that the modifications used in Case 
E3 are sensitive to the grid resolution. This leads us 
to employ the modification (15!))) in Cases E4 and E5. 
Eq. (|56[) is adopted from with some deformation. 
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FIG. 3: The rms of the change in the trace of extrinsic curva- 
ture between consecutive time steps as functions of time with 
different resolutions. 



This modification offers a dissipation effect on Aij which 
fits into our requirement for stability. With this modifica- 
tion, the evolution in Case E4 does not show instability 
until t ~ 950M with AK Ims near round-off error. We 
push the stability in Case E5 by setting the parameter 
used in Eq. (gBJ) to be a = 3/5. Then the AK rms and 
all the other changes drop exponentially until they reach 
round-off error in Case E5. In Case E6, the modification 
(|61~j) instead of Eq. ([56]) is employed. This modification 
gives a little faster convergence. The settings in Case N7 
are the same as in Case E6 except there is not any grid 
symmetry in N7. This case is used to test the effect of the 
grid symmetry on the stability of an evolution with the 
modifications. It shows in Case N7 no sign of instability 
with the relaxation of symmetry. 

Now we would like to look closer at the stable case, i.e., 
Case E6. The results for Case E6 are presented in Fig. [5] 
The upper-left panel shows two different integrations of 
the ADM masses. The (red) dashed line is computed 
from a surface integral (|2"5|) at large separation, while 
the solid line is computed from a volume integral (|29|) 
plus a surface integral over a small sphere enclosing the 
black hole singularity. We choose a radius of i?i = 2M 
for the inner surface and i?2 = 11M for the outer sur- 
face. For R 2 — > oo the two mass integrals should agree 
and should yield the analytic value M of the initial data. 
Our two mass integrals agree to within about 0.05%. The 
lower-left panel in Fig. [5] shows surface and volume inte- 
grations of the angular momentum, similar to the mass 
integrations explained above. The (red) dashed line is 
computed from the outer surface integral ([50")) ; the solid 
line is computed from a combination of volume integral 
([31]) and inner surface integral. For both integrations the 
angular momentum is very close to zero, as it is supposed 
to be. These results indicate that the global quantities 
are consistent with the expected values in single BH and 



are not sensitive to the existence of the excision (as long 
as inside the BH horizon). 

The upper-right panel shows the L2 norms of the 
Hamiltonian constraint % (solid line) and the momen- 
tum constraint M. x (red-dashed line). The lower-right 
panel shows a log plot of the rms of the changes in the 
lapse a (solid line) and the trace of extrinsic curvature 
K (red dashed line) between consecutive time steps. The 
changes in a and K both decrease exponentially until 
they reach round-off error at about t ~ 800 M. We then 
continue to run the numerical simulation for this case till 
the time being over 2000M and there is not any sign of 
instability. 

In Fig. [31 we compare the results of Cases E6's us- 
ing the same modifications but with different resolutions 
to test the convergence of the formulations. From the 
plot we can see that in all the three cases the AK Ims 's 
reach the round-off errors. The cases for the low and the 
medium resolutions are extended over 2000A/ . The case 
for the high resolution is terminated at t ss 1100 A/ due to 
its time-consuming computation. In Case E6-hi we use 
a = 4/5, instead of a = 3/5 in Cases E6-lo and E6, to en- 
hance the stability against higher frequency noise. The 
results of these three cases indicate a nice convergence 
with these modifications. 



B. Binary black hole with punctures 

We use our numerical code AMSS-NCKU, which has 
been employed in several previous studies [HI HH, lil j. 
With this code the standard moving box style mesh re- 
finement is implemented. We used 11 mesh levels in all, 
where 8 levels are fixed and 3 levels are movable. For 
fixed levels we used one box with grids 144 x 144 x 72 
where we have taken the advantage of equatorial sym- 
metry of the system. The outermost physical boundary 
is set to 517.12. For movable levels, two boxes are used. 
And every box has grids 72 x 72 x 36. In time direc- 
tion, the Berger-Oliger numerical scheme is adopted for 
levels higher than 4. Although our code can adopt the 
Kreiss-Oliger numerical dissipation, we disable it in order 
to check the accuracy and stability of our new modifica- 
tion. We set a = 0.1 in Eq. (|4"6l and f(a) = ha in 
Eq. (SI. 

In this subsection both the usual BSSN formalism 
and our modification are tested and compared with each 
other, without the Kreiss-Oliger numerical dissipation for 
both formalisms. Several typical configurations of binary 
black hole listed in the QC sequence of [45j are tested. We 
find that, for QC1 and QC2, both formulations give sta- 
ble and consistent results for whole inspiral and merger 
phases. However, for QC3 and QC6, both formulations 
crash at almost the same time moment. We suspect that 
the failure mainly comes from the numerical error being 
too much. Compared with the tests in the previous sub- 
section, there is one more complicated numerical issue, 
that is, the numerical noise introduced by mesh refine- 
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FIG. 4: Comparison of i — 2, m = 2 mode of ^4 for the usual BSSN formulation and the one with our new modification. The 
outermost boundary locates at 517.12. And the detector locates at R — 50. The (red) dashed line corresponds to the result 
with our new modification and the solid line corresponds to the result with the usual BSSN formulation. The results from the 
two formulations are quite consistent. But for the late ringdown phase, shown in the enlarged subplots, our new modification 
gives a much smoother result. 



merit. And this type of noise causes the numerical insta- 
bility dominantly when the mesh refinement technique is 
employed during the simulations. Compared with the nu- 
merical noise by the discretization of the field equations, 
the numerical noise induced by mesh refinement is much 
larger such that the modification is not able to suppress 
it in time. So in such kind of cases the advantage of our 
modification on stability does not prevail. 

Interestingly we do find our new modification has the 
advantage on accuracy compared with the usual BSSN 
formulation. The advantage is manifest on the gravi- 
tational radiation wave form calculation. And this is- 
sue has been investigated in 0, E(| . As we have men- 
tioned, both the usual BSSN formulation and our new 
modification are stable for simulations of QC1 and QC2. 
For completeness, here we list the initial parameters for 
QC1 and QC2 here. The initial positions of binary 
black hole for these two configurations are (0, ±1.364, 0) 
and (0, ±1.516, 0). The initial linear momentum are 
(T0.286, 0, 0) and (t-0.258, 0, 0). For both configurations 
the two black holes are spineless and identical with punc- 
ture mass parameter 0.463 and 0.47 respectively. We plot 
the resulted waveform in Figf?J The waveform is mea- 
sured by Newman-Penrose quantity \&4. Our calculation 
follows the description in [18[ . The figure shows that the 
two formulations give consistent waveform. And the con- 
sistence convinces us the new modification does work well 
for binary black hole simulations. In particular, the for- 
mulation with our modification suppresses the numerical 
error more effectively and produce a smoother waveform 
than the usual BSSN formulation does during the ring- 
down phase. We should emphasize that the boundary 
effects has no effect on the late-time behavior of wave- 
form. Since the outermost boundary is at 517.12 and the 
waveform is measured at R = 50, therefore the noise in 



the waveform mainly comes from the evolution itself. 

VII. SUMMARY 

We experiment with various modifications of the BSSN 
formulation and study their effects on the stability of nu- 
merical evolution calculations of single black hole. Based 
on the modifications in [13|, we enhance the unimodular 
determinant constraint with Eq. Q36p and the traceless 
extrinsic curvature constraint with Eq. (|37p . We further 
modify the evolution equation for the conformal connec- 
tion functions into Eq. (|49[) by enhancing the linear term 
of r\ With an irreducible decomposition (|3"51) . we re- 
place the component T l g in the connection with the T 1 
and enhance the spatial secondary unimodular determi- 
nant constraint T J \j = to form a new connection T 1 ^ in 
Eq. (|4"Tj). Meanwhile, the field equation off 1 is adapted 
from Eq. (l4"5)) into Eq. (fSTJf for the new connection and 
to ensure the suppression ability of numerical error in 
its linear term. With the irreducible decomposition on 
Tij.fe, the field equation of the conformal metric is modi- 
fied into Eq. PB")) by adding the T-constraint. The field 
equation of the conformal extrinsic curvature is modified 
into Eq. (|61[) by combining the ^-adjustment proposed 
in |14j | and the irreducible decomposition on A^^. 

We found that these modifications on the BSSN for- 
mulation do show their superiority on numerical error 
suppression compared with the earlier work when ap- 
plied to single Kerr-Schild black hole calculations and 
thus increase the accuracy eventually. When applied to 
the binary black hole calculations with the typical ini- 
tial data and without the Kreiss-Oliger dissipation, the 
modified BSSN formulation gives a consistent and more 
accurate result than the conventional method. 
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Among these modifications, Eq. (l46l) seems play a key 
role in stability and thus the whole evolution is quite sen- 
sitive to the value of the parameter a. We suspect that 
this modification, combined with the Ricci curvature in 
the field equation of Aij, could change the characteris- 
tic of T 1 , and thus of the system. However, a further 
analytic/numerical study is needed to have a better un- 
derstanding of the effect of this modification. 

In this work, we only demonstrate the advantage of 
these modifications to stability and accuracy of binary 
black hole simulations without any systematic study on 
the optimal choice of the related parameters. Thus we 
plan to address this problem in more detail in a sepa- 
rate upcoming work to increase the numerical accuracy 



of future binary black hole simulations. 
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